Properties of the Ideal Ginzburg-Landau Vortex Lattice 
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The magnetization curves M{H) for ideal type-II superconductors and the maximum, minimum, 
and saddle point magnetic fields of the vortex lattice are calculated from Ginzburg-Landau theory 
for the entire ranges of applied magnetic fields Hci < H < Hc2 or induction < B < fJ,oHc2 and 
Ginzburg-Landau parameters 2~^^^ < « < 1000. Results for the triangular and square flux-line 
lattices are compared with the results of the circular cell approximation. The exact magnetic field 
B{x,y) and magnetization M{H,k) are compared with often used approximate expressions, some 
of which deviate considerably or have limited validity. Useful limiting expressions and analytical 
interpolation formulas are presented. 

PACS numbers: 74.25.Qt, 74.25.Ha, 74.20.De 



I. INTRODUCTION 

Since Abrikosov'si prediction of the flux-line lattice 
in Type-II superconductors from Ginzburg-Landau (GL) 
theorji^, several approximate formulas for the magneti- 
zation M = B/fiQ — 7? versus the applied magnetic field 
H or average induction B have been publishediiii^i^iSii. 
In these papers and below, the basic situation is con- 
sidered where a macroscopically large, homogeneous and 
isotropic, long superconductor is exposed to a uniform 
parallel field H. In this ideal case demagnetization ef- 
fects, flux-line pinning, and surface effects may be disre- 
garded, and thus the flux-lines are straight lines form- 
ing an ideal periodic lattice. These results are eas- 
ily extended to anisotropic superconductors (where an 
anisotropic effective-mass tensor is introduced into the 
GL theory) by defining an effective GL parameter k that 
depends on the orientation of the flux lines; this trans- 
formation works when H is along a principle symmetry 
axis^i^i^. Generalization to geometries where demagne- 
tization effects occur, are possible by introduction of a 
demagnetizing factor; but this concept works only for 
homogeneous specimens with the shape of an ellipsoid. 
In this case the flux lines in the bulk are still straight 
and form an ideal flux-line lattice (FLL). For other spec- 
imen shapes the FLL is distorted, i.e., the orientation and 
density of the FLL varies spatially and can be calculated 
only numericall y^ ^ i ^ ^ . 

The aim of the present paper is to compare the widely 
used approximate expressions for M{H, k) with the exact 
value obtained numerically and to give useful general an- 
alytic interpolation formulas valid in the entire ranges of 
H and k where the FLL exists, namely, Hd < H < Hc2 
for H, or < B < Bc2 = HoHc2, and 1/V2 
for K, where Hci(T) and Hc2(T) are the lower and upper 
critical fields and k is the GL parameter. Interestingly, 
such general formulas have not been published yet, and 
thus the accuracy of the commonly used expressions is 
not known, probably due to the difficulty of the numer- 
ical solution of the complex-valued GL equations. Early 
numericsii^ used the circular cell method (CCM), which 
approximates the hexagonal unit cell of the triangular 



FLL (or the quadratic unit cell of the square FLL) by a 
circle and the two-dimensional (2D) solution by the ID 
rotationally symmetric solution inside this circular cell; 
both the GL function and the magnetic field are forced 
to have vanishing slope on this circular boundary, as the 
exact solution has on the boundary of the Wigner-Seitz 
cell. This method yields the exact Hd and is expected 
to be best at low inductions B ^ Bc2 where the flux 
lines are well separated. But surprisingly, the circular 
cell approximation gives very good magnetization curves 
at all B (see Fig. 1) and even yields the exact value of 
the upper critical field i?c2- Some more exact results 
of the CCM are listed below. Another methodi uses a 
similar circular symmetric GL order parameter and a lin- 
ear superposition of circular symmetric magnetic fields to 
obtain excellent approximate M{H, k), see also Ref. IT3. 

An in principle exact numerical method— uses pe- 
riodic real trial functions for the squared GL function 
|V'(a;,y)P and magnetic field B{x,y) and minimizes the 
resulting free energy functional with respect to a finite 
number of Fourier coefficients. The same method was 
later applied^® to solve the microscopic BCS-Gor'kov the- 
ory for the properties of the FLL in the entire temper- 
ature interval < T < Tc where Tc is the supercon- 
ducting transition temperature (GL theory, strictly spo- 
ken, applies only close to Tc). Recently this variational 
method was improvedii by keeping the same periodic 
trial functions but now solving the GL equations itera- 
tively; this iteration works much faster and allows to use 
many more Fourier coefficients (many thousands instead 
of only five in Ref. IT5|) . I shall use this 2D iterative pre- 
cision method of Ref. ^3 for the calculation of the FLL 
at i? > 0. At low inductions B <C i?c2 this 2D method 
is supplemented by an iterative circular cell method pre- 
sented in Appendix A. This ID method yields accurate 
values of /ici(k) = Hci/Hc2, which then can be used in 
interpolation formulas. For convenience, I introduce the 
reduced fields b = B/Bc2, h = H/Hc2, m = A//i?c2, such 
that one has m = b — h, hd <h<l, 0<b<l, and 
—hd < m < 0. 

For completeness it should be mentioned that the 
isolated vortesii^ and the FLLi^ have also been com- 
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FIG. 1: Magnetization curves of the triangular FLL, which 
coincide within line thickness with the results for the square 
FLL and for the FLL obtained from the circular cell approx- 
imation, see Fig. 3 for the difference. Shown are h — H/Hc2 
versus b — B/Bc2 (upper left triangle) and — m = —M/Hc2 
versus h (lower right triangle) . One has m = b — h. The lower 
panel shows an enlarged scale. The solid lines show the exact 
numerical result. The dotted lines show the simple interpola- 
tion Eq. (22) good for «: < 5 (upper panel) , and the combined 
low and high field limit Eq. (23) good for k > 1 (lower panel). 



puted from BCS theory (valid at all temperatures) using 
the quasiclassical Eilenberger theory based on energy- 
integrated Green functions. This method was recently 
extended to compute the FLL structure and local den- 
sity of states for s-waveiSiSSiSi, d-wave2iiS^, and chiral p- 
wave2^ superconductors. Very recently the GL methodic 
was generalized phenomenologically to lower tempera- 
tures and to charged vortices24. 



II. TRIANGULAR AND SQUARE FLUX-LINE 
LATTICES AND THE CIRCULAR CELL 
METHOD 

The properties of the FLL within GL theory are cal- 
culated by minimizing the GL free energy of the super- 
conductor with respect to the complex GL function ipir) 
and to the vector potential A(r) of the local magnetic 
field B(r) = V X A. In the usual reduced unita^i^i^i^i^'^''^ 
(length A, magnetic field ^/2Hc, energy density noH^, 
where He is the thermodynamic critical field) the spa- 
tially averaged free energy F of the GL theory referred 
to the Meissner state (■0 = 1, B = 0) reads 



F 



(1) 



Here (...) = (l/V) J cPr . . . means spatial averaging over 
the superconductor of volume V. Introducing the super- 
velocity Q(r) = A — Vip/K and the magnitude /(r) = \ip\ 
of ip{r) = /(r) exp[iif{r)] one may express as a func- 
tional of the real and gauge invariant functions / and 

Q, 



F = 



fQ^ + (V X Q)' 



(2) 



In the presence of vortices Q(r) has to be chosen such 
that V X Q has the appropriate singularities along the 
vortex cores, see e.g. Eq.(B4) in App. B. 

In this paper I consider the ideal periodic FLL in a ho- 
mogeneous (pin-free) large superconductor in a uniform 
magnetic field H along z. In this 2D situation one has 
/ = fix, y),Q^ Qix, y), and B = zB{x, y). Within GL 
theory in reduced units the properties of this ideal FLL 
depend only on two parameters: the GL parameter k and 
the average induction B = {B{x,y)). The equilibrium 
magnetic field H, and the magnetization M = B/hq—H, 
are obtained either from the definition H = dF/dB or, 
more elegantly, from the virial theorem discovered by Do- 
ria, Gubernatis, and Rainer^S^ which in reduced units 
reads 



H = 



{.f-.r + 2B{x,yy 
2B 



(3) 



Some of the properties of the FLL, and all prop- 
erties of the isolated flux line, may be calculated in 
an elegant way by the circular cell approximation2ii^ii4 
as described in App. A. In the circular cell method 
the hexagonal Wigner-Seitz cell around each flux line 
is replaced by a circle with radius R and same area 
ttR^ — ^o/B if each fiux line carries one quantum of 
flux $0 = h/2e = 2.07 • 10"'' Tm^ In reduced units 
one has $o = 27r/K and R/X = R = (2/6«;2)i/2 with 
b = B/Bc2- The boundary conditions on the CCM circle 
r = R are df /dr = dB/dr = 0. I find that the free energy 
of the triangular FLL, Ftr-, and its magnetization, Mtr, 
are reproduced by the CCM with high accuracy in the 
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entire ranges of k and B, 1/^/2 < n < oo and < 6 < 1. 
In particular, the CCM not only yields Hd (in the limit 
i? — » oo) but it also reproduces the exact upper critical 
field Hc2{k), and in the special case k = l/\/2 even the 
exact result H{B) ~ const ^ He ~ Hd = Hc2- These 
somewhat surprising features of this approximation are 
related to the facts that -ffc2, and in the case k = l/\/2 
even the entire curve H{B), are independent of the de- 
tailed arrangement of the flux lines, i.e., they are the 
same for triangular and square or honey-comb FLLs and 
for any other arrangement of single or multiple quanta 
flux lines. Another surprising finding is that the virial 
theorem, Eq. (3), works perfectly in the CCM. Figure 1 
shows the magnetization curves M{H) and the equilib- 
rium field H{B) of the superconductor obtained by the 
CCM for K = 0.85, 1, 1.2, 1.5, 2, 3, 5, 7, 10, and 20. 

In the limit 6^0 the CCM yields the lower criti- 
cal field iJci, which with high accuracy is fitted by the 
formula 



cl 



47rA2 



[ In K 



q;(k) ] 



Hci In K 

q(k) = Qfoo + cxp[— Co — Cl In K — C2(ln k)^] ± e (4) 

with ttoo = 0.49693, Cq = 0.41477, ci = 0.775, Ca = 
0.1303, and e < 0.00076. This expression yields at k = 
I/V2 the correct value hd — 1 and for k ^ 1 it has 
the limit a = 0.49693. A simpler expression for 
yielding an hd with error still less than 1% and with the 
correct limits at k = l/\/2 and k 1, is 



a{K) = 0.5 + 



1 + ln2 

2k - V2 -I- 2 



(4a) 



The CCM in principle cannot yield properties related 
to the different symmetries of the FLL, or to its shear 
modulus, and it cannot give the form factors (Fourier 
coefficients) of the magnetic field B{x, y) that may be 
measured by neutron scattering. These subtle properties 
can be computed by the 2D method presented in Ref. 
and in App. B. This effective numerical method expresses 
the smooth functions /(a;,?/)^ and B{x,y) as 2D Fourier 
series and determines the Fourier coefficients by iteration. 

Figure 2 (top) shows the difference of the free energy 
densities of the triangular {Ftr) and square (Fsq) FLLs. 
This difference is proportional to the shear modulus cqg 
of the triangular FLL (the shear modulus of the unsta- 
ble square FLL is negative within GL theory) by the 
relationii 



C66 



(5) 



Note that this difference is very small, < {Fsq — 
Ptr) / ilJ-oHc) < 0.0018. Even smaller (by ten times) is 
the difference between the free energy densities of the 
CCM (Fee) and of the triangular FLL plotted in Fig. 2 
(bottom). One has < (Fee - Ftr)/inoH^) < 0.00020. 



This result shows that the CCM is an excellent approxi- 
mation for global properties of the FLL. Both differences 
are largest for large k and have a maximum near b ~ 0.3. 
The finding F^q > Ftr means that the triangular FLL is 
stable for all k > 1/V2. Note that for k — 1/V2 one has 
exactly Fgq = Fee = Ftr = for all h. 

Figure 3 (top) shows the difference between the mag- 
netizations Msq of the square FLL and Mtr of the tri- 
angular FLL. Again, this difference is small, 0.0008 < 
-{Msq - Mtr)/Hc2) < 0.00014 and the relative differ- 
ence has the limits -0.018 < {Msq- Mtr) /Mtr < 0.0095. 
Figure 3 (bottom) shows the difference between the mag- 
netization Mcc obtained by the CCM (see Fig. 1) and 
the exact value Mtr of the triangular lattice. Like with 
the free energy, this difference is again smaller by a fac- 
tor of ten than the difference between two lattice sym- 
metries, 0.00016 < -{Mcc - Mtr)/Hc2 < 0.00008 and 
-0.0011 < {Mcc - Mtr) /Mtr < 0.0017. The differences 
vanish exactly at k = 1/V2, and also at k ^ oo, since 
there m = M/Hc2 0. The relative differences (insets 
in Fig. 3) are maximum at k ^ 1. 

The smallness of these differences explains why in 
Fig. 1 the magnetization curves for all three cases Mtr, 
Msq, and Mcc coincide within line thickness. 

Figure 4 shows an example {b = 0.3, k = 1.5) compar- 
ing the spatial functions / and B of the triangular FLL 
and from the CCM. Shown are the cross sections /(a:,0) 
along the nearest neighbor direction x and f{0,y) per- 
pendicular to this, and f{r) from the CCM [a is the vor- 
tex spacing, a^/A^ — An/ {^/Sbn'^)]. It is seen that f{x, 0) 
and f{r), and also B{x,0) and B{r), coincide closely; at 
lower & < 0.3 the difference is smaller than the line thick- 
ness. The solutions for the square FLL deviate more from 
the circular cell solutions. 

The maximum, minimum, and saddle-point fields 
of the triangular FLL, B„iax = ^(0,0), B„„„ = 
B{0,a/^/3), and Bsad = B{a/2,0), depend on b and k. 
Bmax is only slightly above the equilibrium field H, and 
Bsad and Bmin are close to each other and lie some- 
what below the average field B. Bmax and Bmin are 
shown in Fig. 3 of Ref. as functions of b for several 
K = 0.707 ... 5. In Fig. 5 the small differences Bmax — H, 
Bsad — B, and Bsad — Bmin are plotted versus b, in units 
Bc2 and multiplied by a function of k such that the curves 
for all K > l/\/2 collapse at 6 ^ 1. One finds for all k 
near 6=1: 



Bmax H ^ ^^^^^^ K 



-2 - 0.5 



Bc2 



(^2 - 0.069)' 



^(l-&)^ 



B. 



B 



Bc2 
Bsad Bmin 
Br2 



-0.146- 



1 



0.0526 



k2 - 0.069 ' 
k2 _ 0.069 



(6) 
(7) 
(8) 



The factor 0.069 in Eqs. (6)-(8) is 0.5 - 0.5//3a = 0.0688 
where (3 a = 1.1596 is the Abrikosov parameter of the 
triangular FLL. Plots of Bcc{R) — Bmin where Bcc{R) is 
the field value at the boundary of the circular cell in the 
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FIG. 2: Top: The difference of tlie free energy densities of the 
triangular (Ftr) and square (Fsq) FLLs in units HoH^, plot- 
ted versus the reduced induction b = B/Bc2 for k — 0.85, 
1, 1.2, 1.5, 2, 3, 5, 10, and 200. This difference equals 
(2/37r'^) = 0.068 times the shear modulus cge of the trian- 
gular FLL. Bottom: The very small difference between the 
free energy densities of the circular cell method (Fee) and of 
the triangular FLL. Note that the top and bottom plots look 
similar, but the scales of the ordinate differ by a factor of 
about ten. 



CCM, look similar to the plots of Bsad ~ Bmin in Fig. 5 
(lower panel), since the value Bcc{R) lies approximately 
in the middle between Bmin and Bsad, see Fig. 4. Since 
for K 3> 1 and b <C the field in the vortex center 



equals B„ 



2Hci, one has B„ 



H Hci, and 



thus the function plotted in Fig. 5 (upper panel) for & — > 



tends to the limit {bmax ^ h) x I 
0.50), cf. Eq. (4). 

The variance of the magnetic field is 



clK 



^(In K 



a = {[Bix,y)^B]' 



{B{x,vf-B^ 



K#0 



(9) 
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FIG. 3: Top: The difference between the magnetizations 
Msq of the square FLL and Mtr of the triangular FLL in 
units Hc2, plotted versus the reduced induction b = B/Bc2 
for K = 0.85, 1, 1.2, 1.5, 2, 3, 5, 10, and 200. The inset shows 
the relative difference. Bottom: The difference between the 
magnetization M^c obtained by the CCM (see Fig. 1) and the 
exact value Mtr of the triangular lattice. The inset shows the 
relative difference. 



where i?K are the Fourier coefficients of B{x,y) — 
Sk cos Kr and K the vectors of the reciprocal lattice 
of the FLL (App. B). Near 6=1 the Abrikosov solution 
of the linearized GL theor}iiii2& yields for all k values2i 

. = 7.52. 10- S 



B 



A4 (k2 

= 0.172- 



c2 



0.069)2 ' 
1 - b 
- 0.069 



(10) 



The functions 5* and S/{l — b) are plotted in Fig. 6 versus 
Vb for various k. It can be seen that Eq. (10) is a rather 
good approximation for the large range 0.25 < 6 < 1. 
At smaller b the variance a(b) has a maximum and then 
goes to zero again at 6 = 0. 
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FIG. 4: Comparison of the GL functions / and magnetic 
inductions B calculated for the triangular FLL and from the 
circular cell approximation for the example b = 0.3, k = 1.5. 
Shown are the cross sections f{x,0), B{x,0) along the near- 
est neighbor direction x, and f{0,y), B{0,y) along the per- 
pendicular direction y, and /(r), B[r) from the CCM. All 
B are in units _B(0,0) of the triangular FLL. Small devia- 
tions can be seen only close to the cell boundary r = R, 
R/a = 3^/''(27r)"^/^ = 0.525. At lower b the deviations are 
even smaller. 



For small inductions 6 <C 1 and large n one can use 
the London approximation Sk = B/{1 + K'^X'^). For the 
appropriate cutoff at large magnitudes K ~ = k/X 
see Ref. 2B^29. and below. In the range 0.13//t^ < 6 < 1 
the unity in the denominator of Bj^ may be disregarded 
since K^X^ > {in / V3)bK'^ = 7.2556k^ Thus B drops 
out and a becomes independent of b:^'' 



0.00374, S 



(11) 



A4 '■ 

This often used formula corresponds to the upper axis in 
Fig. 6. One can see that this approximation is good only 
for very large k > 70 and only in the range of b near the 
maximum of a. At very small b ^ 0.13/k^ both a{b) 
and S{Vb) drop linearly to zero when b 0. In this 
range the sum in Eq. (9) can be evaluated as an integral, 
yielding 



bn^ $g 



5/2 



(12) 



This approximation is good for k > 5 and very small 6 
(6 < O.OI/k^ for k = 5; 6 < 0.04/^2 for n > 10), see the 
two straight lines in Fig. 6. For k > 5 a good approxi- 
mation, with less than 5% error from 6=1 down to the 
value b ~ 0.25/k^'^ where the maximum of a occurs, is 



1 



B. 



c2 



0.172 ^-[1 + 1.21(1 -%/6)^]. (13) 



This approximation is much better that the interpola- 
tion, Eq. (3) of Ref. 113 



(''max" '"^ ^ ~ 0-069) / (k"^ - 0.5) 



X (1 - b)- 
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FIG. 5: Maximum field Bmax ~ bmaxBc2 minus applied field 
H, saddle-point field Bsad ~ bsadBc2 minus induction B, and 
Bsad minus minimum field Bmin = bminBc.2, for the triangular 
FLL, plotted versus b = B/Bc2 for k = 0.85, 1, 1.2, 1.5, 2, 
3, 5, 7, 10, 20, 50, 100, 200. The solid lines show these small 
differences in units Bc2, multiplied by appropriate functions 
of K to obtain collapse of the curves near b = 1. The dashed 
lines show the same functions multiplied by factors (1 — 6)~^ 
and (1 — b)~^ such that they tend to a finite constant value 
near & = 1, cf. Eqs. (6) - (8). 



III. MAGNETIZATION CURVES 

This section presents analytic expressions which ap- 
proximate the computed magnetization to ~ M/Hc2 — 
b—h (Fig. 1) as a function of the induction b = B/Bc2 or 
of the thermodynamic field h = H/Hc2- We distinguish 
approximations working at high or low inductions. 



A. Approximation for high inductions 

The linearized GL theory yields for 1 — 6 <C 1 
Abrikosov's Bc2 solutioniiii 



1 - b 



(2^2 - 1)/3a + 1 



(14) 



where Pa = {u^DKu^a? = 1 + E™,„exp[if^„^/(47r)] 
(Ref. IT1II30I and App. B) is the Abrikosov parameter, 
Pa = 1.1595953 {Pa = 1.1803406) for the triangular 
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o^'^ X (k^ - 0.069) / B 



X (1 - b)- 




FIG. 6: The magnetic field variance a — {[B{x,y) ~ B]^) of 
the triangular FLL for ft = 0.85, 1, 1.2, 1.5, 2, 3, 5, 7, 10, 20, 
50, 100, 200 plotted in units of Bc2 as ^ ■ {k^ - 0.069)/Bc2 
(solid lines) such that the curves for all k collapse near h = \, 
cf. Eq. (10). The dashed lines show the same functions divided 
by (1 — 6) such that they tend to a finite constant 0.172 at 
6=1. All curves are plotted versus 

\/6= 

to stretch 

them at small h values and show that they go to zero linearly. 
The limits, Eq. (12), for k = 5 and k = 10 are depicted as 
dash-dotted straight lines. The upper frame 0.383 shows the 
approximation (11). 



(square) FLL. The linear magnetization TOyi(6, k) is a 
good approximation in the range 0.5 < 6 < 1, see Fig. 1. 
This suggests the following fit to the exact m: 



- (1 - 6)^exp[/i(&)]gi(K) + ei, 
/i(6) = 2.50u2 - 8.08w + 0.39, w = (1 - 6)° '*\ 

-2-25)(2k2_i)/(2k4). (15) 



m(6, k) 
gi(K) = (1.133+ 1.926/k^ 



with relative error |ei/m| < 0.0013 for h > 0.5 for the 
triangular FLL. Formula (15) is a good approximation 
with relative error < 1% for all k in the large range of 
fields 1/(4^2) + 5 • 10-^ < & < 1, see Fig. 7. 

The same expression (14) fits also the m{b, k) of the 
square FLL, with somewhat larger error if the same func- 
tions /i(6), gii^) are used rather than the optimally fit- 
ted ones. For the difference mtr — Wsq see Fig. 3. 



B. Approximation for "intermediate fields" 

For completeness I mention here also the London 
approximation'^ which was supposed to be good in the 
"intermediate field range" Hd <C -ff ^ Hc2 that exists 
only in superconductors with extremely large k. Within 
London theory the induction is (see App. B) 



K 



cos Kr 



(16) 
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FIG. 7: Lower panel: The exact magnetization M of the 
triangular FLL (solid lines) and the fit, Eq. (15), (solid lines 

with dots) plotted for many k values versus 
to stretch the low field region. Shown is —M normalized to 
its maximum value Hc\ occurring at fe = 0. The fit (15) is 
good for all k and not too small h > 1/{4k^) +0.0005. Upper 
panel: The deviation SM of the fit from the exact M is very 
small when b > 0.5. The dotted lines in the lower panel show 
the old London approximation, Eq. (18). 



where the sum goes over all -FT-vectors with length from 
K = to some cutoff K ^ Inserting this into the 

London free energy density [B(r)^ + A2(V x B)2]/(2/zo) 
and averaging over the superconductor one gets 



F : 



K 



_B2 
2^0 



S$o f d^k 1 
'2^ J 47r2 k^X^ 



(17) 



The integral from k^^^ « (Kw/^f « TT^B/<i>o to k^^^ 

= 2TrBc2/^o equals {4:TtX'^)~'^ \n{-f' /b) where 7' is 
some constant and b = B/Bc2 as above. This yields 



M 



H - 
-M 

Br9 



B 



dF B 
dB ^0 
0.358 



$0 , 7 

79 T 

87rA2^o b 



4k2 



In- 



(18) 



with constant 7 = ^' /e — 0.3575 . . . obtained by our 
fit to the numerical m(6) at k = 200. This old London 
approximation is shown in Fig. 7 as dotted lines. One sees 
that this fit works only at large ft > 20 in the relatively 
small interval 1/(2^^) <b< 0.01, i.e. at very low b (but 
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K = 0.85, 1 .2, 2, 3, 5, 1 0, 20, 50, 200 




FIG. 8: Exact magnetization of the triangular FLL (solid 
lines) and the logarithmic fit, Eq. (19), (solid lines with dots) 
plotted versus h^^^ to stretch the region at small b = BjBci- 
The dotted lines show the London nearest-neighbor approx- 
imation, Eq. (20). The dashed lines show the London ex- 
pression, Eq. (C6), with the sum taken over all shells up to 
Vmax = 100 vortex spacings. Both London approximations 
are good fits at very low b and all k. 



not too low b). It gives to = at 6 = 7 for all k. This fit is 
slightly improved by replacing ln(7/6) by ln(l — 7-I-7/6), 
which gives the correct m = at 6 = 1. 

A much better fit in the spirit of this logarithmic ap- 
proximation is (see Fig. 8) 



m 



1 

4^ 



In 



1 



-/2(6) 



/2(&) = 0.357 -h 2.890 6- 1.58162 



(19) 



This fit is good for k > 3 (error < 3%) and k > 5 (error 
< 1% ) in the large ranges (In/t + 1)/{10k^) < b < 1 
for K = 3 . . . 200. These intervals of validity may also be 
expressed as -M/H^i = -m/hd < 0.8 (0.85) for k < 20 
{k > 50). 



C. Approximation for low inductions 

All the above approximations do not describe the cor- 
rect vertical slope of M{H) aX H = Hd, or zero slope of 
H{B) and miity slope of M(B) = B-H at B = 0. This is 
achieved by the London approximation of pairwise inter- 
acting vortices described in App. C. For very small 6 <S 1 
one may account only for the nearest neighbor shell of six 
vortices in the triangular FLL of spacing a = cX. With 
h{b), Eq. (C8), this yields for -m(6) = h{b) - b: 



m : 



hci 



1 + 



19 



47 



8c 128c2 



a 
A 



(20) 



Formula (20) correctly describes the steep diverging slope 
of m{h) — > 00 or slopes m(6)' — > 1 and h{by — > as 6 — > 
and is valid for < 6 < 2.5/ k'^ for k > 7, Accidentally it 
also fits well m{b) for k < 2 and b < 0.2, see the dotted 
lines in Fig. 8. A smoother fit is obtained by the exact 
London expression (C6) if one or three neighbor shells are 
included in the sum. But taking more terms in the sum 
improves the fit only at large k. Accounting for neighbors 
up to i/ = 100 lattice spacings apart (about 5000 terms) 
one gets good approximation to m and /i for < 6 < 0.01 
(0.02, 0.05) if K > 20 (k = 7, K 2), see Fig.8. In the 
limit 00 the infinite sum (C6) reproduced Eq. (18), 
i.e, the dashed curves in Fig. 8 for k = 50, 200 then will 
straighten and cut the axis M = 0at5 = 7 = 0.358 
(51/2 ^ 0.60 in Fig. 7, b^/^ = 0.71 in Fig. 8). 



D. General interpolation 

All the approximations for m{b) and h{b) known so 
far, including the above formulas, fit either the low or 

high field regions. The formulas (15) and (20) [or better, 
Eq. (C6) with the sum taken over three shells] have a 
small overlap for all k and thus, together, they fit the 
entire range < 5 < 1 [though the good fit of the 1ow-k 
data by the London expression (20) or (C6) is accidental]. 

For practical purposes one may construct interpola- 
tion formulas that approximate the numerically obtained 
magnetization in the entire range < 6 < 1. They should 
satisfy the five conditions 



hiO) = hci, h'{0) = 0, h{l) = 1, h'il) = l-p, 
h"{l) =0, p = m'(l) = [ (2«2 - 1)13^ + 



(21) 



with /ici(«) from Eq. (4). A simple expression that sat- 
isfies all these conditions is 

- to(6, k) =h-b = p{l-b) + {hci - p){l - 6)" (22) 

withr7(fi;) = {l—p)/{hci—p). Formula (22) approximates 

the exact —m(b) well for k < 2 with relative deviation 
|e| < 3%, for K = 3 with -2% < e < 6%, and for k = 5 
with —1% < e < 16%, see the dotted lines in Fig. 1, top. 

For large n, general interpolation formula are more dif- 
ficult to construct because of the nonanalytic limiting 
expression Eq. (20). One may, however, combine the 
miow from Eq. (20) with the TOhigh from Eq. (19) using 
a smooth transition at 6 ~ (2^^)^^, e.g., with weights 



and 



= i + itanh[2.5(26K2 - 1)], or slightly 
better, w=^ + ^ erf [2(26^2 - 1) ], yielding 

m{b, k) = (1 — w) TOiow + w TOhigh • (23) 

This interpolation between expressions (19) and (20) 
works well for < 6 < 1 with relative error |e| < 2% 
for K > 5 and —3.5% < e < 2% for k > 1, see the dotted 
lines in Fig. 1 (bottom). Thus, to in the entire ranges of 
b and k may be approximated by Eq. (22) or Eq. (23). 
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IV. CONCLUSION 

The properties of the ideal periodic flux-hne lattice in 
superconductors are calculated from Ginzburg-Landau 
theory for the entire ranges of GL parameters l/\/2 < 
K < oo and inductions < & = B/Bc2 < 1- The dif- 
ferences between the free energies and magnetizations of 
the triangular and square vortex lattices and the values 
obtained by the circular cell approximation are investi- 
gated in detail. Approximate analytical expressions are 
given for the variance a{b, k) of the periodic induction 
and for the magnetization m{b,n). These limiting and 
interpolation formulas should replace previous approxi- 
mate expressions that have rather limited validity. 

The numerical methods presented in the appendices, 
in principle may be applied also to theories going beyond 
the isotropic GL theory considered here. 



APPENDIX A: ISOLATED VORTEX AND 
CIRCULAR CELL METHOD 

The calculation of the isolated flux line and of the FLL 
within the circular cell method, is a cylindrically symmet- 
ric problem. The free energy depends on the magnitude 
of the GL function /(r) and on the magnetic induction 
B{r) (along z) related to the vector potential A{r) and 
supervelocity Q{r) (along ip) by 



N 



B = 



{ArY {QrY 



Q 



(Al) 



In reduced units V^Hc = noH"^ = A = 1, the free en- 
ergy of a flux line or of the circular cell with radius R 
{■kR^ = <I>o /B) averaged over this cell and referred to the 
Meissner state (/ = 1, i? = 0) reads 



, (/') 



t\2 



+rQ'+B' 



2'Kr dr 



(A2) 



with /' = df /dr. Minimizing the functional (A2) with 
respect to /(r) and Q{r) we obtain the two GL equations, 
which may be written in the form 



f" + K^f = K\2f-f 



-Q^f) + .f'/r, 
B' = fQ = J , 



(A3) 
(A4) 



where j = B' is the current density. In Eq. (A3) a term 
K^f was added on both sides to improve the convergence 
of the iteration below. The boundary conditions are 



f{0) = f'{R)=j{0)=j'{R)=0. 
An appropriate ansatz in terms of Fourier series is 

7r(2m - 1) 



M 



/W=E/GsinGr, G = 



2R 



(A5) 



(A6) 



TV 



A{r) 



aKsmKr+ -B, K=—, 

Tl=l 



(A7) 



X sin. Kr + Kr COS Kr „ 
B{r) = 2_^aK hS, 



N 



Q{r) = aK sin Kr 



1 - rVi?2 



n=l 



N 



Kr cos Kr — {1 + K^r )s\n Kr 



(A8) 
(A9) 
(AlO) 



For the equidistant grid = {i—\)R/Nr, i = 1,2, ... Nr 



one has the orthogonality relation 



sin Gri sin G'n = - Nr 5gg' 



i=l 



(All) 



and similar equations for sin and cos Kri. The GL 
equations (A3) and (A4) thus may be written in the form 
of equations for the Fourier coefiicients fa and qk- 



fa 



1 



G2 + k2 



sin Gri x 



[^\2f-f-Q'f) + f'/r,], 



ax 



N 



1 



K^ + 1 



ax 



— — sin Kri x 



i=l 



(E 



Kr cos Kr — sin Kr 



aK' 



(A12) 



(A13) 



These two equations may be used to obtain the /g and 
aK by iteration, starting with appropriate initial values. 
The iteration becomes more stable and faster if the value 
of the previous iteration step is added with a certain 
weight (1 — c) < 1, e.g., c=0.6, according to the algo- 
rithm: 

/g^(1-c)/g+ cFg{/,Q}, (A14) 
OK ^ (1 - c)aK + cAk{!, 0} , (A15) 

with the symbols FQ{f,Q} and Ak{!,Q\ denoting the 
right hand sides of Eqs. (A12) and (A13), respectively. 
Rapid convergence is achieved by iterating equations 
(A14), (A15) alternately. The equilibrium magnetic field 
H is then obtained from Eq. (3), and the magnetization 
from 



M 



2 

BR ./o 



r - f 



+ B'^-B{rf rdr. (A16) 



At very large k and very small b a large number Nr 
of grid points rj is ne eded to achieve high accuracy, 
Nr ^ R/^ = Rk = \/2/b. In this case the accuracy 
with a limited number of grid points may be improved 
by choosing a nonequidistant grid, e.g., rj = with 
equidistant ui = (i — ^)VR/Nr. To use the orthogonality 
relations one then has to express /, B, and Q as Fourier 
series in the new variable u = and also write the 
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two GL equations in terms of the variable u, using, e.g., 
/'(r) = f'iu)/2u and /"(r) = f"iu)/4u' - f'{u)/4u^. 
This yields 

f"{u) = Au'^^-f + f + Q\f) + f'/u , (A17) 
B'{u) = 2ufQ, (A18) 



and the Fourier series 

M 



/H=E/GsinGw, G 



ni—l 
N 



7r(2m - 1) 
2R 



(A19) 



A{u) ^^axsmKu+^B, A' = ^ , (A20) 



N 



Bin) = 



2 ' R 
2 sin Ku + Ku cos Ku 



2v? 



B, (A21) 



n=l 



n=l 

i^M cos Ku— sin iiTu 



1 — u /i? 

Q(w) = VaKsin/sTii , (A22) 

(A23) 



4m4 



The equations for the new Fourier coefficients are 



fa 



^=l 



G2 + 4^2 

(iu^K^f (1 - /2 - Q2) + /'/^^)] ^ (A24) 



ax 



2 

+ ^ sin Kui 



N 



Ku. COS Ku. - isu^Ku. . (A25) 



For better convergence a term —An fa was added on both 
sides of Eq. (A17) to yield (A24). The corresponding it- 
eration scheme using (A14), (Af 5) needs a smaller weight 
c and more iteration steps, but for large /b it is faster 
than the first scheme since it needs less grid points Nr to 
reach the same accuracy. 



APPENDIX B: PERIODIC VORTEX LATTICE 

The properties of the ideally periodic FLL within GL 
theory may be calculated by minimizing the GL free 
energy of the superconductor, Eq. (2), with respect to 
appropriate periodic trial functions, e.g., Fourier series 
with a large number of terms. For the smooth function 
oj = /^(r) we write the ansatz 



K 



aK(l — cosKr) 



(Bl) 



with r — {x,y), K. — {K^, Ky). In all sums here and be- 
low the term K = is excluded. For vortex positions R = 



R-mn = {mxi+nx2^ ny2) the reciprocal lattice vectors are 
K = Kmn = (27r/5)(mj/2, nxi + 171x2) with 5 = Xiy2 — 
^o/B the unit cell area and m,n ~ 0, ±1, ±2, .... For 
the triangular lattice one has X2 — xi/2, y2 — xi^/3/2, 
and for the square lattice X2 = 0, 1/2 = a^i- For superve- 
locity Q and induction B = V x Q = B{r)z we choose 



B{r) = B + ^&K cosKr, 



(B2) 



K 



z X K 



Q(r) = Q^(r)+^6K-^sinKr. (B3) 



K 



Here ClA{x,y) is the supervelocity of the Abrikosov Bc2 
solution, which satisfies 



V X 



B 



^oYS2{r 

R 



(B4) 



where (52(r) — S{x)6{y) is the 2D delta function. This 
relation shows that is the velocity field of a lattice of 
ideal vortex lines but with zero average rotation. Close 
to each vortex center one has QA(r) ~ z x r'/ {2Kr'^) and 
uj{r) (X r'2 with r' = r — R. In principle Qyi(r) may be 
expressed as a slowly converging Fourier series by inte- 
grating (B4) using divQ = divQyi = as in Ref. 15. But 
it is more convenient to take Qa from the exact relation 



2 K LOA 



(B5) 



where uja^^-, y) is the Abrikosov Bc2 solution given by the 
rapidly converging series (Bl) with coefficients2fli21 



«K = -(-l) 



exp[-X2^„5/(8^)] (B6) 



for general lattice symmetry, and 
— (— 1)''^ exp(— 7rj/2/v^) (y^ ^ m'^ + mn + n'^) for the tri- 
angular lattice. This lua is normalized to {ujA{x,y)) — 1, 
which means that — 1 for any lattice symmetry. 

Another strange property of the Abrikosov solution 
(B6) is that {Vuja/uja)'^ - V^w^/w^ = Att/S = const, 
although both terms diverge at the vortex posi- 
tions; this relation follows from (B4) and (B5) using 
B = ^o/S ^ 27r/{nS). The useful formula (B5) may be 
proven via the complex Bc2 solution ^A{x,y)] it means 
that near Bc2 the third and fourth term in Eq. (2), 
are identical. 

Approximate solutions w(r) and B{y) may be com- 
puted by using a finite number of Fourier coefficients ok 
and 6k and minimizing the free energy k, ok, 6k) 
with respect to these coefficientsi^. However, a much 
faster and more accurate solution methodic is to iter- 
ate the two GL equations 5F/5lo ~ and SF/SQ, = 
written in appropriate form. Namely, the iteration 
is stable and converges rapidly if one isolates a term 
(^V^ -I- const) (w, Q) on the l.h.s. and puts the remaining 
terms to the r.h.s. as a kind of "inhomogeneity" of such 
London-like equations, e.g., 

(-V^ + 2k2)cj = 2k,^ {2u; - iv^ - LjQ^ - g), {B7) 
{-W^+u;)Qb = - ojQa - {co ~ c^)Qb , (B8) 
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with the abbreviations g{r) — (Woj)'^ /{Ak^uj), Qh = 
Q - Qa> V X Qb = B{r) - B, and u! = (uj) = 
X^K^K- Equations (B7), (B8) introduced some "pen- 
etration depths" (2k2)-i/2 ^ ^/y2 and Q-'^/^ = A/w^/^ 
(in real units), which stabilize the convergence of the it- 
eration. Acting on the Fourier series oj (Bl) and Q;, (B3) 
the Laplacian yields a factor —K^; this facilitates the 
inversion of (B7) and (B8). Using the orthonormality 

2 (cos Kr cos K'r) = ^kk' (B9) 

(for K ^ 0) one obtains from (Bl), (B2) ok = 
— 2(a;(r) cosKr) and &k = 2(i?(r) cosKr). The conver- 
gence of the iteration is considerably improved by adding 
a third equation which minimizes F, Eq. (2), with respect 
to the amplitude of uj, i.e., dF/dui — 0. This step gives 
the largest decrease of F. The resulting three iteration 
equations for the parameters ck and 6k then readii 



flK := 



4k2((w2 + ojQ'^ _2uj + g) cos Kr) 



if 2 + 2k2 

-2([(w-cj)B(r) -hp]cosKr) 



3K 



(BIO) 
(Bll) 
(B12) 



with p = (Vo; X Q)z = Q^duj/dy — Qyduj/dx and g — 
{\/Lof/{4.K^w) = (V/)Vk2 as above. 

The solutions ijj(r), B(r), and Q(r) are obtained by 
starting, e.g., with ok = (1 — b) and = and 
then iterating the three equations (BIO), (Bll), (B12) by 
turns until the coefficients do not change any more. Af- 
ter typically 25 such triple steps, the solution stays con- 
stant to all 15 digits and the GL equations are exactly 
satisfied. Since all terms in (BIO) - (B12) are smooth 
periodic functions of r, high accuracy is achieved by us- 
ing a regular spatial 2D grid, e.g., Xi = [i — l/2)xi/Nx 
(i = 1 ... TV,) and y, = (j - l/2)y2/(2iV,) (j = 1 • . . Ny, 
2Ny w Nxy2/xi) with constant weights xi/N^ and 
y2/{2Ny). These N = N^Ny = 100 to 5000 grid points 
fill the rectangular basic area 0<a:<a:i, 0<y< 2/2/2, 
which is valid for any unit cell with the shape of a par- 
allelogram. Spatial averaging (...) then just means sum- 
ming N terms and dividing by N . 

Best accuracy is achieved by considering all Km„ 
vectors within a half circle |Kmn| < Kmax, with 
^max ~ 20N/S chosen such that the number of the 
Kmn is slightly less than the number N of grid points. 
The high precision of this method may be checked with 
the identity B(x,y) / Bc2 = 1 — ^{x,y), which is valid 
at K = I/V2 for all b. This relation is confirmed with 
an error < 10^^. The equilibrium field H or reversible 
magnetization M = B — H is computed from Doria's 
virial theorem, Eq. (3). 

APPENDIX C: LONDON THEORY 



where 6{x,y) is the 2D delta function. The solution for 
the magnetic field of one isolated vortex at R = is 



B,(r) = ($o/2^A2)Xo(r/A). 
The modified Bessel function 

(Pk coskr 



i^o(./A)=/-- 



(C2) 



(C3) 



has the derivative ii'o(a;)' = —Ki{x) with the limits 
Kq{x^I) « -Inx, Ki{x^l) « 1/x, and for a;>l(S 

r ^ - J, I 9 225 A 



For a periodic FLL one obtains the Fourier series i?(x, y), 
Eq. (16), which may also be written as a sum over isolated 
vortex fields, B{x,y) = ^p^i3^(r — R). Similarly, the 
free energy of the FLL may be written as a sum of vortex 
self energies ($o-f^ci per unit length) plus a double sum 
over all interactions between two vortices. The average 
energy density F, Eq.(17), then reads 



BHr 



R, 



(C5) 



For the triangular vortex lattice we write R/X — vc with 
c = a/A = (47r/A/3)^/^(&K2)~^/^ (a = vortex spacing) 
and v"^ = w? + mn + = 1, 3, 4, 7, 7, 9, ... . Taking 
the derivative H = dF/dB one obtains for h = H/Hc2 
with hci = Hci/Hc2- 



h = hci + ^V] IXoii^c) + ^Ki 



{vc) 



(C6) 



Here the sum is over v — \^ \/3, 2, . . . , i.e. the number of 
six flux lines per shell is already accounted for. Equation 
(C6) is still exact. It works for 5 <C 1 (i.e. for nonover- 
lapping vortex cores) and for k > 1.4 (i.e. when the long- 
range interaction of vortices is purely magnetioiiiS) . 
With the expansions (C4) one obtains for x = vc~^ 1: 



19 47 

IH ^ 

8a; 128x2 



(C7) 



At very small &, namely for c = a/A 3> 1, the sum may 
be restricted to the nearest neighbor shell, i.e. to the first 
term, v = \, yielding 



The modified London equation for a lattice of straight 
vortex lines at regular positions R — R„i,i (App. B) is 



(1 



AV^) 



B(x,y) 



<i>o2^<5(r 

R 



Rri 



(CI) 



hc\ -I- 



2k2 



1 1^ _ 47 
^ 8c ^ 128c2 



(C8) 
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